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Abstract 

A  study  is  performed  to  find  the  optimal  operating  conditions  of  hydrogen  polymer  electrolyte  fuel  cells  using  an  efficient  optimization 
approach  based  on  validated  multi-resolution  fuel  cell  simulation  tool  developed  in  house.  Through  the  design  of  experiment  method,  a  set 
of  designed  simulation  runs  were  carried  out  using  the  fuel  cell  simulation  tool.  Based  on  the  simulation  results,  an  analytic  metamodel 
was  then  constructed  using  the  radial  basis  function  approach.  A  feasible  sequential  quadratic  programming  scheme  was  then  employed  to 
optimize  the  metamodel  to  achieve  the  global  optimal  solutions.  To  illustrate  the  optimization  approach,  four  control  parameters  including  cell 
temperature,  cathode  stoichiometry,  cathode  pressure,  and  cathode  relative  humidity  were  considered.  The  optimization  objective  is  defined 
as  the  maximization  of  the  overall  efficiency  of  the  fuel  cell  system  under  ideal  or  realistic  system  assumptions.  The  study  shows  that  different 
optimal  solutions  exist  for  different  system  assumptions,  as  well  as  different  current  loading  levels,  classified  into  small,  medium,  and  large 
current  densities.  The  approach  adopted  in  this  study  is  generic  and  can  be  readily  applied  to  a  larger  number  of  control  parameters  and  further 
to  the  fuel  cell  design  optimizations. 

©  2005  Elsevier  B  .V.  All  rights  reserved. 
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1.  Introduction 

Numerous  components  and  control  parameters  are 
involved  for  a  hydrogen  polymer  electrolyte  fuel  cell  (PEFC) 
system  to  operate  at  the  optimal  level.  Previous  studies  [1,2] 
on  the  balance  of  plant  of  fuel  cell  systems  focused  on  the 
optimization  of  systems  consisting  of  fuel  cell  stacks,  com¬ 
pressors,  humidifiers,  and  cooling  units.  These  optimization 
studies  used  phenomenon  fuel  cell  models  that  were  obtained 
through  thermodynamic  analysis  and/or  empirical  data.  How¬ 
ever,  due  to  the  existence  of  various  irreversible  processes 
inside  fuel  cells  and  the  limited  availability  of  empirical 
data,  the  phenomenon  models  often  yield  incorrect  predic¬ 
tions  when  the  empirical  data  is  not  available,  especially 
under  the  operating  conditions  with  various  gas  humidity 
levels  and  stoichiometries.  Recent  experimental  study  [3]  on 
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hydrogen  PEFCs  showed  that  the  performance  of  fuel  cells 
could  be  greatly  affected  by  the  anode  and  cathode  humid¬ 
ity  levels,  stoichiometries,  cell  temperatures,  and  different 
combinations  of  these  conditions.  Meanwhile,  different  fuel 
cell  designs  and  materials  used  in  the  fuel  cell  can  also  have 
significant  effects  on  the  fuel  cell  performance  [3].  Since 
phenomenon  models  are  difficult  to  cover  various  fuel  cell 
operating  conditions  and  fuel  cell  designs,  a  model  that  cap¬ 
tures  the  physics  associated  with  each  of  the  processes  in  a 
fuel  cell  has  to  be  employed  for  more  accurate  predictions  of 
PEFC  operations. 

Theoretically,  the  existing  full  three-dimensional  (3D) 
approaches  [4,5]  can  be  employed  to  study  the  details  of 
PEFCs,  but  simulations  of  a  full  size  fuel  cell  and  fuel  cell 
stacks  are  prohibitive  due  to  the  unbearable  computational 
cost.  Further,  for  the  purpose  of  design  and  operation  opti¬ 
mization  of  fuel  cells,  the  simulation  turn  around  time  has 
to  be  dramatically  reduced  to  satisfy  the  industrial  needs. 
On  the  other  hand,  some  simplified  models  such  as  the 
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Nomenclature 

a 

active  area  density  (m2  m-3) 

A 

area  (m2) 

4/o 

active  catalyst  area  x  exchange  current  density 
(A  m-3) 

c 

species  mass  concentration  (kmolm-3) 

D 

diffusivity  (m2  s-1) 

E 

cell  potential  (V) 

F 

Faraday  constant  (96,487  C  mol-1) 

G 

Gibbs  free  energy  per  unit  time  (J  s-1) 

g 

molar  Gibbs  free  energy  (Jkmol-1) 

h 

molar  water  latent  heat  (J  kmol-1) 

i 

local  current  density  (A  m-2) 

I 

average  current  density  (A  m-2) 

K 

permeability  of  porous  GDF  (m2) 

L 

characteristic  length  of  the  agglomerate  (m) 

m 

molar  fraction 

m 

gas  channel  mass  flux  rate  (kgm-2  s-1) 

Ns 

number  of  species 

P 

pressure  (Pa  or  atm) 

R 

gas  constant  (Jkmol-1  K-1) 

t 

time  (s) 

T 

temperature  (K) 

u 

velocity  (ms-1) 

V 

species  diffusion  velocity  (ms-1)  or  cell  volt¬ 
age  (V) 

w 

power  (Js-1) 

y,  z 

Cartesian  coordinates  (m) 

Y 

species  mass  fraction 

Greek  letters 

a 

anode  transfer  coefficient 

p 

cathodic  transfer  coefficient 

8 

thickness  of  surrounding  Nation  layer  (m) 

£ 

porosity  of  gas  diffusion  layer 

0 

membrane  proton  potential  (V) 

Y 

ratio  of  air  specific  heats  (1.4) 

T] 

active  overpotential  (V)  or  efficiency 

[1 

viscosity  (kgm-1  s-1) 

P 

density  (kgm-3) 

a 

effective  ionic  conductivity  (S  m-1) 

i 

stoichiometry 

Subscripts 

a 

anode 

avg 

average 

c 

cathode 

h2 

hydrogen 

hum 

humidification 

k 

the  kth  species 

m 

membrane 

02 

oxygen 

RH 

relative  humidity 

sat 

saturation 

sys 

system 

wall 

gas  channel  wall 

empirical  models  [6],  one-dimensional  (ID)  models  [7,8], 
two-dimensional  (2D)  models  [9-11],  and  quasi-3D  model 
[12],  cannot  be  used  for  systematic  study  of  various  flow  and 
coolant  channel  design  and  optimizations  due  to  the  neglect 
of  various  channel  patterns. 

In  this  paper,  the  multi-resolution  approach  developed 
by  Wu  and  Liu  [13]  was  employed  for  fuel  cell  simulation 
and  optimization  of  operating  conditions  for  PEFCs.  In  the 
multi-resolution  approach,  a  3D  model  is  employed  for  the 
membrane,  in  which  the  transport  includes  both  convections 
and  diffusions  in  all  directions.  The  catalyst  layer  is  modeled 
using  a  ID  +  2D  model,  in  which  at  each  location  of  the  fuel 
cell  plate,  the  governing  equations  are  integrated  only  in  the 
direction  (ID)  perpendicular  to  the  fuel  cell  plate  (2D).  The 
gas  diffusion  layer  is  modeled  as  a  3D  flow  model  due  to  the 
dominance  of  diffusion  process.  Modeling  the  flow  channels 
presents  the  most  challenges,  because  a  full  size,  high  power, 
single  cell  can  have  very  long  and  complex  flow  channels.  In 
the  multi-resolution  approach,  the  flow  channels  are  repre¬ 
sented  by  quasi- ID  models  with  conservation  laws  satisfied 
at  the  interface  between  the  flow  channel  and  the  gas  diffusion 
layer. 

Although  the  multi-resolution  simulation  approach 
reduces  the  simulation  time  to  a  small  fraction  of  that  of 
the  full  3D  approach,  the  computational  time  is  still  too  high 
to  directly  combine  a  optimization  procedure  with  the  simu¬ 
lation  program,  because  doing  so  would  require  hundreds  or 
even  thousands  of  sequentially  performed  simulations.  To  this 
end,  the  metamodeling  approach  has  to  be  adopted  to  reduce 
the  number  of  simulations  and  make  optimization  feasible. 

In  the  metamodeling  approach,  an  approximate  function 
(metamodel)  is  created  to  replace  the  simulation  program  in 
optimization  to  obtain  an  output  response  for  a  given  input. 
Since  creating  the  metamodels  only  requires  a  limited  num¬ 
ber  of  predefined  simulations  and  evaluating  the  metamodels 
is  very  efficient,  the  total  cost  for  obtaining  the  optimal 
designs  is  very  small.  Polynomial  regression  has  been  the 
most  common  metamodeling  method  for  years,  but  they  are 
only  appropriate  for  linear  and  quadratic  responses.  Recent 
studies  showed  that  radial  basis  functions  (RBF)  were  more 
accurate  in  creating  models  for  both  low-  and  high-order 
non-linear  responses  [14-17].  In  this  study,  due  to  the  high 
non-linearity  of  the  fuel  cell  responses,  the  RBF  models  are 
used  in  our  optimization  of  fuel  cell  operations. 

Since  this  study  focuses  on  the  systematic  approach 
to  achieving  the  optimal  operating  conditions  for  PEFCs 
through  multi-resolution  simulation,  the  feasibility  of  the 
approach  will  be  demonstrated  through  the  optimization  of 
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operating  conditions  of  a  single  PEFC  cell  with  control 
parameters  including  cell  temperature,  cathode  stoichiom¬ 
etry,  cathode  gas  pressure,  and  cathode  relative  humidity. 
The  optimization  with  the  objectives  of  achieving  the  highest 
fuel  cell  efficiencies  is  performed  using  the  multi-disciplinary 
optimization  software  developed  by  Fang  and  Horstemeyer 
[18]. 

In  the  remaining  portion  of  this  paper,  the  theoretical 
background  of  the  multi-resolution  simulation  is  first  given. 
A  brief  review  on  the  RBF  is  then  presented  followed  by 
the  design  optimization  formulations.  The  simulation  setup 
and  the  optimization  objectives  are  then  explained  in  detail. 
Finally,  the  validation  of  the  simulation  models  and  the  anal¬ 
ysis  of  the  optimization  are  performed,  followed  by  the  con¬ 
clusion. 


Table  1 


Cathode  catalyst  layer  parameter  values 


Symbol 

Parameter 

Values 

4/0 

Specific  exchange  current  density 
(Am-3) 

5  x  104 

f 

Cathodic  transfer  coefficient 

2 

(7 

Effective  ionic  conductivities  of 

Nation  (S  m-1) 

7 

8/a 

Thickness  of  surrounding  Nation 
layer  divided  by  active  agglomerate 
area  density  (m2) 

6.2  x  10“12 

L 

Characteristic  length  of  the 
agglomerate  (m) 

3  x  10“6 

r\eff  1 1 

^02,Nafioii 

Effective  oxygen  diffusivity  in 

Nation  times  Henry’s  law  constant 
(m2  s-1) 

4  x  io-15 

2.  Multi-resolution  fuel  cell  simulation 


oxygen  at  the  cathode  catalyst  layer.  The  local  current  density 
of  the  anode  catalyst  layer  za  is  modeled  by  the  Butler- Volmer 
equation  as  [4] 


In  this  section,  the  agglomerate  oxygen  reduction  model 
employed  for  the  cathode  catalyst  layer,  the  interface  model 
for  the  anode  catalyst  layer,  the  porous  media  model  for  both 
gas  diffusion  layers  (GDF)  and  the  membrane  layer,  and  the 
quasi- ID  model  for  the  flow  channels  are  presented  along 
with  the  related  interface  conditions. 

2.7.  Agglomerate  catalyst  layer  model 

The  electrochemical  reactions  in  the  cathode  catalyst  layer 
can  be  symbolically  expressed  as 

02  +  4H+  +  4e“  2H20.  (1) 


H  =  i 


ref 

0,a 


CH2 

z^ref 

CH2 


Vh2 


exp 


—  exp 


acF  V 
RT  ria)_ 


where  ig  is  the  anode  reference  exchange  current  den¬ 
sity;  Cg  the  anode  hydrogen  reference  concentration;  yn? 
the  hydrogen  concentration  parameter;  aa  the  anode  anodic 
charge  transfer  coefficient;  ac  is  the  anode  cathodic  charge 
transfer  coefficient.  The  parameters  used  in  Eq.  (3)  are  listed 
in  Table  2. 


2.3.  Gas  diffusion  layer  model 


A  first  order  oxygen  reduction  mechanism  is  employed  in 
association  with  the  agglomerate  model,  in  which  the  catalyst 
layer  is  treated  as  the  collection  of  agglomerates  consisting  of 
carbon  particle  support  and  platinum  particles  surrounded  by 
a  thin  layer  of  Nation.  In  contrast  to  the  macro-homogeneous 
model,  which  treats  the  catalyst  layer  as  a  homogenous  phase, 
the  agglomerate  model  is  capable  of  predicting  a  voltage 
drop-off  at  large  current  densities  by  integrating  the  mass  dif¬ 
fusion  transfer  resistance  [13,19].  The  parameters  and  their 
values  used  in  this  model  are  given  in  Table  1 . 

2.2.  Anode  catalyst  layer  interface  model 

The  electrochemical  reactions  in  the  anode  catalyst  layer 
are  given  by 

H2  -*  2H+  +  2e~ .  (2) 

Hydrogen  is  transported  to  the  catalyst  layer  through  the 
GDF,  and  discharges  electrons  by  the  electrochemical  reac¬ 
tion  described  in  Eq.  (2).  The  electrons  move  through  the 
external  circuit  to  provide  useful  current,  while  the  protons 
transport  through  the  polymer  electrolyte  membrane  to  the 
cathode  catalyst  layer  and  produce  water  by  the  reduction  of 


The  GDF  model  consists  of  the  classic  Darcy’s  law  for  the 
velocity  field  and  the  gas  phase  species  transport.  The  species 
transport  equations  are  derived  from  those  of  the  mass  con¬ 
servation  laws.  The  governing  equations  of  the  GDF  model 
are  given  as  follows: 


dp 

£  +  V  •  ( pu )  —  0, 

ot 

(4) 

u  =  — 

00 

< 

1 

(5) 

dPYk 

£ - 

dt 

+  V  •  (puYk)  =  -V  •  {pYkVk), 

(6) 

Table  2 

Anode  catalyst  layer  parameter  values 

Symbol 

Parameter 

Values 

,-ref 

*0,a 

Anode  reference  exchange  current  density 

6  x  103 

(Am-2) 

Anode  hydrogen  reference  concentration 
(kmolm-3) 

1.2 

z^ref 

ch2 

Kh2 

Hydrogen  concentration  parameter 

0.5 

aa 

Anode  anodic  charge  transfer  coefficient 

1 

ac 

Anode  cathodic  charge  transfer  coefficient 

1 
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where  s,  p ,  u ,  K ,  fi,p,  g ,  and  F&  are  the  porosity,  gas  phase  den¬ 
sity,  super-facial  velocity,  permeability,  effective  viscosity, 
pressure,  gravity,  and  the  kth  species  mass  fraction,  respec¬ 
tively.  The  diffusion  velocity  Vk  of  the  kth  species  is  modeled 
by  Fick’s  law  as 

YkVk  =  -DfVYk.  (7) 

where  Dff  is  the  effective  kth  species  diffusivity  and  is  cal¬ 
culated  using  the  Bruggeman  relation  as 

Df  =  eL5Dk.  (8) 


component,  non-reactive  gas  mixture  can  be  written  as 


d 
d  t 


Q  dv  +  F2A2  —  F\A\ 


(13) 


where  A 1  and  A2  are  inlet  and  outlet  areas  and  F\  and  F2  are 
fluxes  through  the  flow  areas.  The  vectors  Q ,  F,  and  S  are 
given  by 


~  pYk  ~ 

pYku  -  pDkVYk 

Q  = 

P 

F  = 

pu 

pu 

pu2  +  p  —  fiWu 

Since  temperature  variations  in  a  single  cell  are  often  neg¬ 
ligible,  the  isothermal  condition  is  assumed  and  the  energy 
equation  is  reduced  to  a  constant  temperature  condition.  The 
transport  coefficients,  such  as  the  viscosity  /z  and  the  gas  dif¬ 
fusivity  Dk  are  computed  using  CHEMKIN  developed  at  the 
Sandia  National  Laboratories  [20]. 

2.4.  Membrane  model 


rhk  d A 


^wall 


(k=l,Ns 


fhk  d A 


1) 


(14) 


where  Awan  is  the  wall  area.  Eq.  (14)  takes  into  account  the 
effects  of  the  species  mass  transport  and  the  channel  viscous 
friction. 


The  water  transfer  inside  membrane  is  modeled  to  account 
for  the  convection  driven  by  pressure  gradient,  electro- 
osmotic  drag  of  water  from  the  anode  to  the  cathode,  and 
diffusion  from  high  to  low  concentration  of  water  content. 
Due  to  the  electro-neutrality  assumption  and  the  homoge¬ 
neous  distribution  of  charge  sites,  constant  proton  concen¬ 
tration  c™+  in  the  membrane  is  assumed.  The  momentum 
equation  takes  the  form  of  generalized  Darcy  relation.  The 
potential  equation  is  derived  from  Ohm’s  Law.  The  system 
equations  for  the  membrane  assuming  isothermal  condition 
are  given  as  follows: 


drm 

acu2o 

dt 


+  V  •  AfH2o  —  0’ 


Km 

Mm  =  ~£ - (V/7  -  pg), 

dm 


l 


F 


m  .  -_cm 


V0m  = - 1 - Cu+U, 


(9) 

(10) 

(11) 


Net  water  molar  flux  in  the  membrane  is  given  by  the 
Nernst-Planck  equation  along  with  the  Nernst-Einstein  rela¬ 
tionship  as 

^H20  =  ^H20VcH20  +  cH20*L  (12) 

where  is  the  drag  coefficient  and  is  the  water  self¬ 
diffusion  coefficient. 


2.5.  Quasi- ID  channel  model 

The  gas  channel  model  is  developed  from  the  mass  and 
momentum  conservation  laws  with  variables  changing  only 
in  the  flow  direction.  The  governing  equations  for  the  multi¬ 


2.6.  Boundary  and  interface  conditions 

The  inlet  velocities,  temperatures,  and  species  mass  frac¬ 
tions  of  the  anode  and  cathode  need  to  be  specified.  At  the 
outlets  of  the  gas  channel,  only  the  backpressure  is  speci¬ 
fied  with  an  extrapolated  boundary  condition  employed  for 
other  variables.  Non-slip  boundary  condition  is  specified  at 
all  external  surfaces. 

For  the  bipolar  solid  wall,  a  non-slip  boundary  condition 
is  employed,  and  no  species  mass  transport  is  allowed.  The 
detailed  formulation  of  the  quasi- 1 D  boundary  conditions  can 
be  found  in  Ref.  [13]. 

At  the  interface  between  the  GDL  and  the  gas  channel,  the 
GDL  provides  the  mass  flux  including  both  mass  convection 
and  diffusion,  and  the  gas  channel  provides  the  GDL  with 
the  values  of  dependent  variables  including  the  pressure  and 
species  concentrations.  The  viscous  friction  is  calculated  the 
same  way  as  the  solid  wall,  because  the  velocity  at  the  inter¬ 
face  is  almost  zero  and  much  smaller  than  the  mean  channel 
velocity. 

The  interface  conditions  between  the  GDL  and  the  catalyst 
layer  are  implemented  as  the  following.  For  the  GDL,  the  kth 
species  mass  flu xfk  including  both  the  mass  convection  and 
diffusion  through  the  interface  is  given  as 

(pu  Yk  -  Df V  Yk)  I  GDL  =  fk  ■  (15) 

At  the  anode  side,  the  hydrogen  consumption  flux  based  on 
the  local  current  density  is  given  as 

M\ 

M  =  -ffw  (16) 

where  ia  is  the  anode  catalyst  local  current  density.  Water  flux 
at  the  interface  is  given  as  Eq.  (12).  At  the  cathode,  oxygen 
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flux  due  to  reaction  is  given  as 


3.  Metamodeling  and  optimization 


(17) 


where  the  local  current  density  ic  is  calculated  through  the 
cathode  agglomerate  model.  Besides  the  water  transfer  due 
to  electro-osmotic  drag,  water  created  due  to  electrochemical 
reaction  at  the  catalyst  layer  is  given  as 


(18) 


The  total  water  flux  at  cathode  catalyst  layer  interface  is  then 
given  as 


/h2o  =  Sh2o  +  A^h2o-  (19) 

It  is  assumed  that  the  water  in  the  GDL  is  in  equilibrium  with 
that  in  the  membrane  layer  at  the  interface,  which  implies 
that  the  water  activity  is  same  for  both  membrane  and  GDL 
at  the  interface. 

The  membrane  potential  equation  is  solved  by  obtaining 
the  local  current  density  at  the  cathode  side  together  with  the 
potential  value  at  the  anode  side.  The  cathode  current  density 
is  calculated  through  the  cathode  agglomerate  model,  while 
the  potential  at  the  anode  side  is  set  as  the  anode  activation 
overpotential.  In  many  previous  studies,  the  activation  over¬ 
potentials  of  anode  and  cathode  are  specified  as  constants. 
But  in  this  paper,  they  are  computed  and  updated  iteratively, 
which  results  in  more  accurate  prediction  of  the  local  current 
density. 

The  fuel  cell  reversible  open  circuit  voltage  is  calculated 
as  [21] 

£re y  =  1.23  -  0.9  X  10“3(T  -  298.15) 

+  |^(ln  /+,  +  I  In  P02)  (20) 

where  Pu2  and  Po2  are  the  partial  pressure  of  hydrogen  at 
the  anode  inlet  and  the  partial  pressure  of  oxygen  at  the  cath¬ 
ode  inlet,  respectively.  After  the  initialization  of  the  cathode 
activation  overpotential,  an  initial  cathode  current  density 
is  obtained.  Further,  the  anode  activation  overpotential  is 
updated  according  to  Eq.  (3).  The  membrane  potential  equa¬ 
tion  is  solved  to  obtain  the  membrane  ohmic  overpotential 
77 ohm  and  finally  the  cathode  overpotential  is  updated  as 


0c  —  ^cell  ^rev  Osl  Oohmi  (71) 

where  £ceii  is  the  specified  cell  potential. 

This  algorithm  is  capable  of  predicting  the  average  cell 
current  density  based  on  the  specified  cell  potential  or  cell 
potential  based  on  the  specified  average  cell  current  density. 
Either  the  flow  rates  or  the  stoichiometry  can  be  specified  for 
the  simulation.  The  data  flowchart  among  the  related  mod¬ 
ules  for  predicting  the  average  cell  current  density  based  on 
specified  cell  potential  and  stoichiometry  is  shown  in  Fig.  1. 


Optimization  is  a  reverse  engineering  process  to  find 
the  input  parameters  corresponding  to  the  optimal  output 
responses.  This  process  is  iterative  and  typically  involves  a 
large  number  of  calculations  to  obtain  output  responses  based 
on  given  inputs.  For  fuel  cell  applications,  the  input-output 
relationships  are  not  available  in  explicit  form  and  must 
be  obtained  through  simulations  that  are  computationally 
expensive.  To  make  optimization  feasible  and  efficient,  the 
metamodeling  approach  has  to  be  adopted  to  represent 
the  input-output  relationships  with  approximation  functions 
(metamodels)  that  are  in  explicit  form  and  less  expensive  to 
evaluate.  Optimization  is  then  performed  on  metamodels  to 
find  the  optimal  solutions.  Since  creating  metamodels  only 
requires  a  limited  number  of  predefined  simulations  (sam¬ 
pling  points),  the  total  computational  cost  is  very  small. 

3.1.  Metamodeling  with  radial  basis  functions 

The  RBF  was  originally  developed  by  Hardy  [22]  to  fit 
large  and  non-linear  data  sets  and  was  shown  to  produce 
good  approximations.  However,  recent  studies  showed  that 
the  RBF  was  also  appropriate  for  both  low-  and  high-order 
non-linear  responses  with  limited  samplings  [14-17]. 

An  augmented  RBF  model  has  the  form  of 

n  p 

f'(x)  =  \\x  -  Xi  ||)  +  ffcjgjix),  (22) 

i= 1  J=1 

where  n  is  the  number  of  sampling  points,  x  a  vector  of  design 
variables  (input  parameters),  Xi  a  vector  of  design  variables 
at  the  /th  sampling  point,  |  |x  —  Xi  \  |  the  Euclidean  norm,  0  a 
basis  function,  A/  the  coefficient  to  be  determined  for  the  ith 
basis  function,  g(x)  a  linear  polynomial  function,  p  the  total 
number  of  terms  in  the  polynomial,  and  cj  is  the  coefficient  to 
be  determined.  In  this  study,  the  Wu’s  compactly  supported 
function  03  o  [23]  is  used  in  creating  all  the  response  functions 
and  is  given  by 

03jO(?)  =  (1  -  07( 5  +  35 1  +  10 It2  +  147/3  +  101 14 

+35 15  +  5 16)  (23) 

where  t  is  the  normalized  Euclidean  norm 
(t=\\x  —  Xi  \  |/max(|  |x£  —  x/|  |),  k=l,  2,  . . .,  n).  Eq.  (22) 
is  underdetermined,  because  there  are  more  parameters  to  be 
solved  than  the  number  of  equations  created  with  available 
data  points.  Therefore,  the  orthogonality  condition  is  further 
imposed  on  coefficients  A  so  that 

n 

y ^Xigj(Xj)  =  0,  for  j  =  1,  2, . . . ,  p.  (24) 

i= 1 

By  replacing  x  and /(x)  in  Eq.  (22)  with  the  n  vectors  of  design 
variables  and  corresponding  function  values  and  combining 
with  Eq.  (24),  we  obtain  ( n  +  p)  equations  in  the  matrix  format 
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Fig.  1.  Data  flowchart  among  related  modules  for  the  multi-resolution  simulation. 


as 


where  Aij  =  </>(\\Xi-Xj\\)  (/=  1,  2,  . . n,j=  1,  2,  . . n ), 
Gij  0"  1  ?  2,  •••?  Hi  j  1  ?  2,  . . . ,  pf  X  1 1  X2,  •  •  •  i 

Xn)T ,  C  =  [CU  C2i  •  • Cp]T,  and /=  \f(x i),/(x2),  •  •  -if(xn)]T. 
Solving  Eq.  (25)  gives  the  coefficients  X  and  c  for  the  RBF 
function  in  Eq.  (22). 

3.2.  Optimization  formulation  and  solution 

The  optimization  problem  of  this  study  is  to  achieve  opti¬ 
mal  operating  conditions  in  terms  of  PEFC  efficiencies  under 
small,  medium,  and  large  electrical  currents,  respectively. 
Four  control  parameters  (design  variables)  are  selected  in 


this  study;  they  are  the  cell  temperature,  cathode  stoichiome¬ 
try,  cathode  gas  pressure,  and  cathode  relative  humidity.  The 
optimization  problem  consists  of  three  tasks  each  of  which 
can  be  formulated  as 

Findx  which  maximizes /(x)  subject  to  the  constraints 

x\<Xi<xf  i  =  1,4  (26) 

where  fix)  is  the  efficiency  function  obtained  through  meta¬ 
modeling,  and  x\  and  are  the  lower  and  upper  bounds  of 
the  ith  design  variable,  respectively. 

The  gradient-based  feasible  sequential  quadratic  program¬ 
ming  (FSQP)  method  [24]  is  used  for  this  optimization 
problem.  Since  the  RBF  models  for  the  objectives  are  typ¬ 
ically  highly  non-linear  functions,  the  FSQP  method  can 
easily  be  trapped  into  local  optima  with  a  single  start¬ 
ing  point.  To  avoid  this  problem,  a  number  of  random 
starting  points  are  used  in  solving  each  of  the  optimiza- 
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Fig.  2.  The  grid  of  the  three-channel  serpentine  gas  flow  channel  employed 
for  both  the  anode  and  cathode  sides. 


tion  tasks,  and  the  best  solution  is  chosen  as  the  final 
optimum. 


4.  Simulation  setup 

The  simulations  were  carried  out  for  a  typical  25  cm2 
PEFC  with  three  channel  serpentine  flow  patterns  for  both 
the  cathode  and  anode  sides  with  co-flow  arrangement.  Fig.  2 
shows  the  gas  channel  grid,  generated  using  the  paramet¬ 
ric  grid  generation  code  developed  in  house.  The  channel 
height,  width,  rib  width,  and  number  of  inlet  and  channel 
turns  are  used  as  parameters  for  the  specific  grid  generations. 
The  channel  height,  width,  and  rib  width  are  set  as  1 ,  1 ,  and 
1.2  mm,  respectively.  The  grid  size  for  the  cathode  GDF,  the 
membrane  layer  and  the  anode  GDF  are  all  10x51  x51.  The 
overall  active  area  of  the  fuel  cell  is  50  mm  x  50  mm.  The 
physical  parameters  of  GDF  and  membrane  layer  (Nafion 
115)  used  in  this  simulation  are  listed  in  Table  3. 

To  ease  the  analysis,  four  parameters,  including  the  cath¬ 
ode  pressure,  relative  humidity,  stoichiometry,  and  the  cell 


Table  3 


Physical  parameters  of  GDL  and  membrane  layer 


Symbol 

Parameter 

Values 

£ 

GDL  porosity 

0.5 

K 

GDL  permeability  (m2) 

1.0  X  10"12 

8 

GDL  thickness  (mm) 

0.3302 

£m 

Membrane  porosity 

0.5 

Km 

Membrane  permeability  (m2) 

1.8  x  10"18 

Pm 

Membrane  effective  viscosity  (kgm-1  s-1) 

8.91 

Pdry 

Dry  membrane  density  (kgm-3) 

2000 

Mm 

Membrane  molecular  mass  (kg  kmol- 1 ) 

1100 

nm 

H+ 

Proton  concentration  (kmolm~3) 

1.2 

<5m 

Membrane  thickness  (mm) 

0.1778 

Table  4 


Operating  parameter  ranges  and  values 


Symbol 

Parameter 

Range  or  value 

T 

Cell  temperature  (°C) 

50-90 

/c 

Cathode  stoichiometry 

1.1-5 

G 

Anode  stoichiometry  (fixed) 

1.5 

Pc 

Cathode  pressure  (atm) 

1-5 

Pa 

Anode  pressure  (atm)  (fixed) 

1 

RHC 

Cathode  relative  humidity  (%) 

10-100 

RHa 

Anode  relative  humidity  (%)  (fixed) 

100 

temperature,  were  chosen  to  adjust  the  fuel  cell  operations, 
while  the  anode  pressure,  relative  humidity  and  stoichiometry 
were  fixed  at  1 .0  atm,  100%  and  1 .5,  respectively.  The  ranges 
of  each  parameter  shown  in  Table  4  were  selected  based  on 
the  previous  experience  of  fuel  cell  operations.  Using  the 
advanced  design  of  experiment  method,  instead  of  thousands 
of  simulation  runs,  only  75  simulation  runs  were  conducted 
in  parallel  for  the  optimization  of  operating  conditions. 

Many  possible  optimization  objectives  exist  for  different 
applications.  Though  the  current  approach  is  generic  with 
various  objectives,  for  the  purpose  of  simplifying  the  anal¬ 
ysis  and  discussion,  the  optimization  objective  of  this  study 
is  chosen  as  the  maximization  of  the  fuel  cell  efficiencies 
under  the  average  current  density  loadings  of  0.15,  0.45,  and 
0.75  A  cm-2  to  represent  small,  medium,  and  large  current 
density  loadings,  respectively. 

The  efficiency  ijsys  of  the  system,  consisting  of  a  single 
cell,  humidifier,  and  air  compressor,  is  defined  as 


Osys  — 


A/ayg  Vcell 


G  +  lOhum  T  ^compressor 


(27) 


where  A,  7avg,  and  Vceii  are  the  fuel  cell  active  area,  aver¬ 
age  current  density  and  cell  voltage,  respectively.  The  maxi¬ 
mum  work  done  per  unit  time  through  consumption  of  equal 
amount  of  hydrogen  under  reversible  conditions  is  the  Gibbs 
free  energy  per  unit  time  G  expressed  as: 


A/avg<?H2 

2F[£a  +  (fa  —  1)?7r] 


(28) 


where  gH2>  fa,  and  OR  are  the  hydrogen  molar  Gibbs  free 
energy,  anode  stoichiometry  and  hydrogen  recycle  efficiency, 
ibhum  and  ^compressor  are  the  power  done  to  humidify  the 
anode  and  cathode  inlet  streams,  and  to  compress  the  air  from 
ambient  pressure  to  cathode  inlet  pressure  [25],  respectively: 


A/avg/ZL 

Z?hum^ 


^c^satFHc  ^a/^satFHa 

4Q?c  —  /?satRHc)-T^  2(/?a  —  i^satRHa) 

(29) 


_  A/avgCpTambient 

^compressor  —  ~T~Z  rjT 

4  F ^compressor 


(— ) 

\  P ambient  / 


y/y- 1 


-  1 


(30) 
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where  £c,  /zl,  /?sat>  -^air?  77hum?  dnd  ^compressor  are  the  cathode 
stoichiometry,  molar  latent  heat  of  water,  water  saturation 
pressure  at  the  cell  temperature,  molar  fraction  of  oxygen 
in  dry  air,  humidifier,  and  isentropic  compressor  efficiency, 
respectively.  In  addition,  RHC  and  RHa  are  the  cathode  and 
anode  relative  humidity,  respectively.  Here,  zhcompreSsor  is 
computed  assuming  the  mechanical  efficiency  of  the  com¬ 
pressor  to  be  unit  according  to  Larminie  and  Dicks  [25]. 
Since  the  anode  side  conditions  are  fixed  in  current  study, 
the  hydrogen  compression  work  required  is  not  considered 
and  also  hydrogen  is  assumed  to  be  fully  recycled.  In  addi¬ 
tion,  it  is  assumed  that  the  only  useful  work  done  by  the  fuel 
cell  is  electric  power  and  no  heat  co-generation  is  considered 
for  the  low  temperature  operation  of  fuel  cells. 

5.  Results  and  discussion 

The  single  phase  multi-resolution  fuel  cell  simulation 
framework  developed  in  house  has  been  validated  against 
the  experimental  data  under  different  operating  conditions 
when  flooding  is  not  a  serious  factor.  Similar  to  the  approach 


used  by  Springer  et  al.  [7],  the  liquid  water  is  treated  as 
super  gas.  As  shown  in  Fig.  3,  the  numerical  predictions  and 
the  experimental  measurements  for  different  cathode  relative 
humidities  are  in  good  agreement.  The  details  of  the  valida¬ 
tion  and  accuracy  analysis  of  the  simulations  can  be  found  in 
Ref.  [26]. 

5.7.  Fuel  cell  performance  under  different  operating 
conditions 

The  fuel  cell  might  be  operated  at  different  range  of  aver¬ 
age  current  densities  for  different  applications.  For  example, 
if  the  fuel  cell  efficiency  and  the  stable  performance  of  fuel 
cells  are  required,  fuel  cells  might  be  loaded  at  lower  current 
density  range.  If  the  high  power  output  is  required,  fuel  cells 
might  be  operated  under  medium-range  current  density  load¬ 
ing.  However,  if  high  current  output  is  required,  fuel  cells  will 
tend  to  be  operated  close  to  the  limiting  current,  constrained 
of  course  by  the  stability  of  the  system. 

Under  different  current  loadings,  the  effects  of  each  oper¬ 
ating  parameter  on  the  fuel  cell  performance  are  different. 
For  example,  from  Fig.  4,  showing  the  effect  of  the  cathode 
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Fig.  3.  Comparison  of  polarization  curves  and  membrane  ohmic  overpotential  for  different  cathode  inlet  air  RH,  while  maintaining  anode  at  80%  relative 
humidity.  Cathode  inlet  air  RH  (a)  100%;  (b)  70%;  (c)  50%;  (d)  30%. 
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Fig.  4.  Comparison  of  polarization  curves  for  different  cathode  inlet  air  RH,  Fig.  6.  Comparison  of  polarization  curves  for  different  operating  pressures 

while  maintaining  a  fully  humidified  anode  inlet  condition.  and  cathode  inlet  air  relative  humidity. 


relative  humidity  on  fuel  cell  performance,  high  cathode  rel¬ 
ative  humidity  tends  to  reduce  the  fuel  cell  performance  at 
large  current  density  when  excessive  water  is  generated  and 
cannot  be  removed  quickly.  On  the  other  hand,  at  medium 
current  density,  when  the  sufficiency  of  the  oxygen  and  the 
membrane  hydration  status  generally  determine  the  fuel  cell 
performance,  the  high  relative  humidity  can  either  adversely 
affect  the  fuel  cell  performance  due  to  the  decrease  in  oxy¬ 
gen  supply  or  favorably  due  to  the  increase  in  membrane 
conductivity.  At  low  current  density,  when  ohmic  poten¬ 
tial  is  negligible,  the  high  relative  humidity  often  tends  to 
adversely  affect  the  fuel  cell  performance  due  to  low  oxy¬ 
gen  partial  pressure,  though  at  very  small  magnitude.  In 
addition,  different  combinations  of  the  operating  conditions 
also  play  an  important  role  in  the  fuel  cell  performance.  If 
the  cathode  oxygen  stoichiometry  is  increased  from  2  to 
10,  the  effect  of  cathode  relative  humidity  on  cell  perfor¬ 
mance  changes  as  shown  in  Fig.  5,  since  the  sufficiency 
of  oxygen  is  not  a  factor  under  high  stoichiometry  and 


Fig.  5.  Comparison  of  polarization  curves  for  different  cathode  inlet  air 
RH,  increasing  the  cathode  stoichiometry  to  10,  while  maintaining  a  fully 
humidified  anode  inlet  condition. 


only  the  membrane  hydration  status  determines  the  fuel  cell 
performance. 

Though  parameters  such  as  the  relative  humidity  may 
have  different  effects  on  fuel  cell  performance  under  different 
conditions,  the  increase  of  cathode  pressure  increases  the  fuel 
cell  performance  in  general  as  shown  in  Fig.  6.  However,  the 
increase  of  cathode  pressure  comes  at  the  cost  of  work  done 
by  the  compressor.  In  reality,  the  favorable  effect  of  high 
cathode  pressure  can  also  be  compromised  by  the  compres¬ 
sor  noise  level  and  gas  sealing  requirements,  which  are  not 
considered  in  the  current  study  for  ease  of  analysis,  though. 

5.2.  Optimization  of  fuel  cell  operating  conditions 

Since  detailed  simulations  of  auxiliary  components  such 
as  humidifier  and  air  compressor  are  not  considered,  as  they 
are  not  the  focus  of  the  current  study,  their  efficiencies  are 
assumed.  To  evaluate  the  effects  of  auxiliary  component 
efficiency  on  the  overall  system  efficiency,  an  ideal  system, 
in  which  the  work  done  by  all  the  auxiliary  components  are 
neglected,  is  first  assumed.  Then,  to  simplify  the  analysis  but 
without  losing  generality,  the  humidifier  and  air  compressor 
efficiencies  are  assumed  both  at  50%.  Table  5  shows  the 
optimization  results  of  the  ideal  case.  As  expected,  the 

Table  5 


Optimization  results  under  ideal  condition 


Parameter 

Current  density 

0.15  A  cm-2 

0.45  A  cm  2 

0.75  A  cm  2 

Cell  temperature  (°C) 

51.1 

61.8 

90 

Cathode  stoichiometry 

3.88 

5 

5 

Cathode  pressure  (atm) 

3.28 

5 

5 

Cathode  relative  humidity 
(%) 

10.6 

10 

75.2 

Efficiency  (true  value) 

0.618 

0.492 

0.393 

Efficiency  (prediction) 

0.619 

0.495 

0.412 

Error  (%) 

0.31 

0.67 

4.68 

J.  Wu  et  al.  /  Journal  of  Power  Sources  156  (2006)  388-399 


397 


optimal  efficiency  drops  as  the  current  loading  increases  due 
to  the  increase  in  total  overpotential  that  consists  of  acti¬ 
vation,  ohmic,  and  concentration  contributions.  In  general, 
the  activation  overpotential  only  increases  slightly  with  the 
change  of  current  density  while  the  ohmic  overpotential 
increases  linearly  when  the  membrane  hydration  status 
does  not  vary  much,  and  the  concentration  overpotential 
will  increase  dramatically  near  and  at  the  limiting  current 
density.  Since  no  compressor  cost  is  involved,  the  fuel  cell 
operation  optimization  is  achieved  near  or  at  the  higher 
bound  of  the  specification  ranges  of  the  cathode  pressure 
and  stoichiometry  for  all  three  levels  of  current  density 
loadings. 

The  optimal  cell  temperature  lies  in  the  lower  bound  of 
the  specification  range  for  small  and  medium  current  density 
loadings  due  to  the  fact  that  fuel  cell  open  circuit  voltage  and 
thus  the  cell  voltage  increase  as  the  cell  temperature  decreases 
from  Eq.  (20).  However,  at  large  current  densities,  the  fuel  cell 
voltage  is  more  influenced  by  the  transport  properties  inside 
the  fuel  cells  and  as  a  result,  the  optimal  cell  temperature  for 
the  large  current  density  loading  is  at  the  higher  bound  of 
the  specification  range,  in  which  the  flooding  effect  is  more 
likely  to  be  reduced  due  to  high  water  saturation  pressure  at 
high  temperature.  In  addition,  low  cathode  relative  humidity 
optimizes  the  fuel  cell  efficiency  at  small  and  medium  cur¬ 
rent  densities  while  relative  high  relative  humidity  optimizes 
the  fuel  cell  efficiency  at  large  current  density  shown  also  in 
Table  5.  The  result  is  actually  the  complex  interplay  of  all  the 
control  parameters.  At  small  and  medium  current  density,  the 
increase  in  partial  pressure  of  oxygen  greatly  increases  the 
fuel  cell  performance  through  reducing  the  activation  over¬ 
potential;  however,  at  large  current  density,  since  the  optimal 
cell  temperature  is  shown  to  be  at  the  high  end  of  the  spec¬ 
ification  rage,  the  membrane  hydration  status  and  thus  the 
proton  conductivity  plays  a  more  important  role:  the  increase 
in  relative  humidity  tend  to  hydrate  the  membrane  and  thus 
reduce  the  ohmic  overpotential. 

To  validate  the  optimization  results,  the  simulations  with 
the  predicted  optimal  parameter  values  were  carried  out 
and  the  true  system  efficiencies  were  then  obtained  for  the 
three  different  current  loadings,  respectively.  The  comparison 
between  the  true  efficiency  and  predicted  optimal  efficiency 
is  also  shown  in  Table  5  for  the  ideal  case.  As  we  can  see 
that  the  maximum  prediction  error  is  only  4.68%  for  the 
large  current  loading  case.  Generally  speaking,  the  prediction 
error  depends  on  the  non-linearity  of  the  response  functions 
and  can  be  improved  by  increasing  the  number  of  sampling 
points. 

Table  6  shows  the  optimization  result  of  the  case  in  which 
the  air  compressor  and  humidifier  efficiencies  are  both  set  at 
50%.  Understandably,  the  system  efficiencies  for  all  the  load¬ 
ings  have  dropped  compared  to  the  ideal  case,  and  especially 
for  the  large  loading  case,  of  which  the  system  efficiency 
has  dropped  around  35%.  It  can  also  be  seen  that  the  opti¬ 
mal  values  for  the  four  control  parameters  have  changed 
greatly.  Due  to  the  consideration  of  the  work  done  by  the 


Table  6 


Optimization  results  under  realistic  condition 


Parameter 

Current  density 

0.15  A  cm"2 

0.45  A  cm  2 

0.75  A  cm-2 

Cell  temperature  (°C) 

50.2 

55.5 

70.4 

Cathode  stoichiometry 

1.34 

1.53 

1.85 

Cathode  pressure  (atm) 

1.69 

2.46 

2.94 

Cathode  relative  humidity 

(%) 

15.9 

10.4 

12.6 

Efficiency  (true  value) 

0.523 

0.374 

0.258 

Efficiency  (prediction) 

0.522 

0.374 

0.265 

Error  (%) 

-0.15 

0.09 

2.67 

air  compressor  and  humidifier,  the  optimal  values  for  the 
cathode  stoichiometry,  cathode  pressure,  and  cathode  relative 
humidity  tend  to  be  at  the  lower  bound  of  the  specification 
ranges.  The  variation  of  optimal  values  for  the  cell  tem- 


Fig.  7.  System  efficiency  iso-surface  plots  under  large  current  density  load¬ 
ing  with  the  cell  temperature  fixed  at  the  predicted  optimal  values  for  the 
ideal  and  realistic  system  cases,  respectively,  (a)  Ideal  systems  and  (b)  real¬ 
istic  systems. 
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perature  with  the  current  loading  levels  still  maintains  the 
same  trend  as  the  ideal  case  though  the  optimal  cell  tem¬ 
perature  value  at  large  current  density  loading  drops  from 
the  higher  bound  to  the  middle  of  the  specification  range. 
Similar  to  the  ideal  case,  the  comparison  between  the  true 
value  and  the  prediction  of  the  system  efficiency  shows  that 
the  optimization  result  is  accurate  with  maximum  error  at 
only  2.67%  and  the  number  of  simulation  runs  should  be 
sufficient. 

For  the  purpose  of  illustrating  the  difference  between  the 
ideal  and  realistic  response  functions,  Fig.  7  shows  the  three 
level  iso- surfaces  of  the  ideal  and  realistic  response  functions 
under  large  current  loading  with  the  cell  temperature  fixed 
at  the  predicted  optimal  values,  respectively.  Essentially,  the 
iso-surface  plots,  with  x-,  y-,  z- axis  being  the  cathode  stoi¬ 
chiometry,  cathode  pressure,  and  relative  humidity,  respec¬ 
tively,  exhibit  the  path  toward  the  optimal  solutions  for  each 
case.  It  can  also  be  seen  that  the  objective  function  (system 
efficiency)  for  the  ideal  case  under  large  current  density  load¬ 
ing  is  much  more  irregular  than  the  corresponding  realistic 
case,  which  also  explains  the  larger  prediction  error  for  the 
ideal  case. 


6.  Conclusion 

An  efficient  and  systematic  approach  to  achieving  the  opti¬ 
mal  operating  conditions  for  PEFCs  has  been  developed.  The 
validated  multi-resolution  fuel  cell  simulation  tool  developed 
in  house  has  been  utilized  to  carry  out  a  specified  number 
of  simulations  designed  by  using  the  design  of  experiment 
method.  Due  to  the  high  irregularity  of  the  model  describing 
the  overall  fuel  cell  system  efficiency  as  a  function  of  the  four 
control  parameters  including  cell  temperature,  cathode  stoi¬ 
chiometry,  cathode  pressure,  and  relative  humidity,  the  radial 
basis  function  approach  for  constructing  the  approximation 
function  or  metamodel  is  employed  instead  of  the  polyno¬ 
mial  regression  approach  that  is  only  suitable  for  regular 
models.  In  order  to  achieve  the  global  optimization,  the  fea¬ 
sible  sequential  quadratic  programming  method  is  modified 
by  starting  from  a  set  of  randomly  picked  initial  points.  The 
overall  system  efficiency  defined  under  either  ideal  or  realis¬ 
tic  system  conditions  is  used  as  the  objective  function.  Three 
typical  levels  of  current  loadings  0. 15, 0.45,  and  0.75  A  cm-2 
are  considered  to  represent  small,  medium  and  large  current 
loadings. 

The  fact  that  the  predicted  optimization  results  are  in  good 
agreement  with  the  true  values  demonstrates  the  accuracy  and 
validity  of  the  approach.  Different  optimal  solutions  exist 
for  different  system  assumptions  as  well  as  different  current 
loadings.  Under  the  ideal  system  assumptions  in  which  the 
work  done  by  the  auxiliary  components  are  neglected,  the 
highest  efficiencies  under  different  current  loadings  tend  to 
be  achieved  near  or  at  the  higher  bound  of  the  specification 
ranges  of  the  cathode  stoichiometry  and  pressure.  However, 
under  the  realistic  conditions  in  which  the  work  done  by  the 


humidifier  and  the  air  compressor  is  considered,  the  opti¬ 
mal  choice  for  the  four  control  parameters  to  achieve  the 
highest  system  efficiency  changes  and  the  maximum  system 
efficiency  also  drops. 

The  approach  developed  herein  is  neither  restricted  to 
only  four  control  parameters  nor  to  the  fuel  cell  operation 
optimizations.  The  approach  will  be  readily  applied  to  the 
fuel  cell  design  optimizations  in  the  future.  When  the  number 
of  control  parameter  increases,  especially  for  the  design 
optimization  applications,  the  design  of  experiment  approach 
can  ensure  that  the  number  of  simulation  runs  needed  to 
stay  within  reach.  However,  an  approximate  estimation  of 
required  number  of  simulation  runs  as  a  function  of  number 
of  control  parameters  is  needed  to  achieve  certain  accuracy 
level. 
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